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Abstract 

The effects of erosion, avalanching and random precipitation are captured 
in a simple stochastic partial differential equation for modelling the evolution 
of river networks. Our model leads to a self-organized structured landscape 
and to abstraction and piracy of the smaller tributaries as the evolution pro- 
ceeds. An algebraic distribution of the average basin areas and a power law 

relationship between the drainage basin area and the river length are found. 
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A fractal river network is a striking example of self-organized criticality. The physics 
of river network evolution arises from an interplay of the structured landscape governing 
the water flow with the erosional effects of the water feeding back into further sculpting of 
the landscape. Extensive studies of the fractal characteristics of real river networks have 
been carried out. |I]-[7| Hack has studied the relationship between the length of a river / 
and the area of a drainage basin s. s is a measure of the total area of the land covered by 
the principal stream and its tributaries that feed into the network. Hack's measurements 
indicate that for basin areas s ranging over almost five decades (up to 375 square miles), 
s ~ r with the exponent 1/0 ~ 0.57 . Other measurements of the distribution of drainage 
basin areas suggest a power law scaling of the form P(s) ~ s~ T with r = 1.45 ± 0.03 . [p]] 

Most of the models of river networks fall into two categories. The first is restricted 
to reproducing the statistical properties of networks ||. More recently, models for the 
evolution of river networks have been developed. Based on careful studies of river data [|TJ , 
Rinaldo and coworkers |J have suggested that the effects of local optimal rules equivalent to 
critical erosion parameters lead to statistical characteristics for the networks similar to global 



constraints of minimum energy dissipation. Leheny and Nagel [|T(J introduced a lattice model 



that incorporated erosion and showed a competition in growth between neighboring river 



basins relevant for late stages of evolution. Kramer and Marder [|ll| have also constructed 
a lattice model that allows for the elucidation of the scaling properties of the large scale 
features of river networks. In addition, they have proposed coupled differential equations 
for two scalar fields, the height of the soil and the depth of the water flowing over the soil. 
An analysis of these equations has led to an understanding of the shape and stability of 
individual river channels. 

Our principal goal is to introduce and numerically study a simple stochastic partial 
differential equation for the evolution of the landscape of a river network. Our model is a field 
theory for the soil height h(x,t) and takes into account the effects of random precipitation, 
erosion and the avalanching of soil. An initially smooth landscape evolves into a nontrivial 
spatially self-organized state in which Hack's law and the algebraic distribution of drainage 



basin areas are obtained. 



The evolution equation may be written in the compact form: 



d_ 
at 



h(x, t)=V- h(x, t) - kV 4 h(x, t) + r}(x, t) 



(1) 



where V- = D\&1 + D 2 (\ Vh \)d% and x = (x,y). The coefficient D 2 {\ V/i |) is defined as: 



and the random noise rj(x,t) = sr(x,t) is Gaussianly distributed with zero average and 
correlation 



The equation describes the temporal evolution of a two dimensional landscape with 
periodic boundary conditions in the x direction and a dominant water flow in the y direction 
due to an incline. The equation allows for ordinary diffusion in the y direction and a 
stabilizing V 4 /i term working as an ultraviolet regulator. In the x direction, a diffusion term 
is operational as long as | Vh | > M. These diffusional processes mimic the avalanching effect 
discussed by Leheny and Nagel: when neighboring sites have large differences in heights, 
avalanching of the soil occurs leading to a smoothing effect. The erosional processes are 
captured by the negative diffusion coefficient D 2 when | V/i | < M. Such a term accentuates 
the height difference and captures the physics of erosion - the water flowing in the shallower 
parts of the landscape erode the soil further, increasing the height difference. In the simple 
version of our model, the erosional processes are blind - the erosion is uncorrelated with 
the water flow. Also, our analysis has been carried out in the limit of a large underlying y 
slope so that the water flow is directed downhill and ordinary diffusion is operative in the y 
direction. Large basins are known to have smaller overall slopes so that a less directed water 
path would be more appropriate. The noise 7] mimics the erosional effects of random 
precipitation. It should be noted that this noise is an essential ingredient of the dynamics: 






(3) 
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noisy initial conditions would evolve into a flat landscape under only the deterministic part 
of equation ([[]) . 

Our equation is a generalization of equations developed in other contexts. In the limit 
that D\ = D2 > 0, the equation reduces to the Edwards- Wilkinson equation for linear 
growth processes ||12||. When Di and D 2 are both positive but different, the equation is akin 



to Barenblatt's equation JT^| describing the pressure in an elasto-plastic porous medium that 



allows for the contraction and expansion of the medium due to the flow. The novel ingredient 
in our equation is the possibility of a negative diffusion coefficient. As in sandpile models 
of self-organized criticality |14]], smoothing diffusional processes kick in when the gradient 



exceeds a critical value. While the sandpile models are driven by the random addition of 
sand, here the instability is caused by the erosion and the noise. Our model evolves in a 
self-organized manner to a noisy boundary separating stability and instability, as do the 
sandpile models. In the sandpile case, there is a well-defined interval between the addition 
of grains during which the relaxation takes place. In contrast, the dynamics is continuously 
turned on in our model. 

In order to discretize eq. ([]]) it is convenient to work with dimensionless variables t,x,h 
defined as: 



* = (4) 



h = e -=-h (6) 

where Ax is the lattice spacing of the spatial mesh we will be using. In the following we 
shall omit the tildas. The discretized version of the equation is: 

h(x, t + At)= h(x, t) + -^L(A y h(x, t)+D[\ Vh \]A x h(x, t)) (7) 
At A 2 h(x,t) + VAt r(x,t) 



Ax A 



Here A x and A y are discrete second derivatives in the x and y direction and A 2 is the discrete 
square laplacian defined so that the Fourier transform of A 2 h does not contain anisotropic 
terms to the leading order, and 

| 1 if I V/i |> M 

D[\ Vh\] = { (8) 

De|<0 if | V/i |< M 

Once the landscape reaches a self-organized state, we obtain a measure of s by adding a 
test drop of precipitation on each site and following its downhill descent along the steepest 
path p|JT0||. A drop at site (x,y) has three possible destination sites - (x — l,y+ 1), {x,y + 1) 
or (x + 1, y + 1). The site with the smallest h is selected. On a given site, s is defined as 
the number of drops collected on that site during the downhill evolution. 

If we define P(s, L) as the density of sites with a total flow s on a landscape of linear 
size L, one may combine Hack's law and the algebraic scaling of the river basin areas into a 



generalized scaling form p5[ : 



P(s,L)=s-TF(s/L*') (9) 

Assuming that rivers are fractals with length I ~ L dp , the characteristic basin area s c ~ 
jj>' i<t>'/ d F implying Hack's law with = 0' /dp. In the present case we find dp = 1, i.e. 
the rivers are self-affine, and thus = 0'. 

From the definition of s, it follows that its average J2 S P{ s i L)s = (L + l)/2 which 
combined with @ leads to the scaling law = (2 — r) _1 . |18| Further, for r > 1, the 
normalizability of P(s,L) in the L —>■ oo limit imposes a constraint that \im x ^o F(x) is 
universal and equal to (r— 1). Figure |] shows a plot of P(s, L) for different sizes L = 10,30,50 
and 100 at time t =10. These results were obtained in the regime where At/ Ax 2 ~ 10~ 3 , 
well below the linear stability limit At /Ax 2 ~ 0.25. The figures are then collapsed into 
a single scaling plot (Figure 0) with the choices of r = 4/3 and = (2 — r) _1 = 3/2. 
p~9|] Within the error estimates of our analysis, these exponents are the same as obtained 
in Scheidegger's static model H for river networks |16||. They are also within the range 



observed for real rivers in the small basin limit. It should be mentioned that the above 



scaling form is expected to be valid only in the large s and large L regime. Long time 
runs on small systems are suggestive of a non-Scheidegger universality class at large times. 
However, a quantitative study of this regime appears to be beyond our present capability 
due to computational limitations. 

We have also monitored the temporal evolution of the roughness of the landscape. The 
roughness W is defined as: 

w(t) = (^j:<(hm-h(t)r >) 1 / 2 

X 

where h(t) = Ylsh(x,t)/L. The average was performed over 300 samples corresponding to 
different realizations of the noise. We find an intermediate regime with an algebraic growth 
W(t) ~ W with an exponent (3 = 0.21 ± 0.02 (Figure This exponent value is different 
from the exactly solvable Edwards- Wilkinson model |12| where (3 = 0, i.e. W(t) ~ Vint. At 
longer times saturation due to finite size is observed. In this context, it is interesting to note 
the recent work of Czirok et. al. ]2(J on a geomorphological micromodel of mountains. The 
evolution of river networks is shown in Figure |], where only sites with s > 30 are shown, 
while the corresponding elevation profile is shown in Figure [5]. The model captures some of 
the key ingredients of dock's theory of the evolution of river networks || - small tributaries 
are eliminated as the main rivers swell in size. The figures also show the effects of piracy in 
which a larger aggressive stream captures the flow of a neighboring stream. Such effects had 
also been shown to occur by Leheny and Nagel |T(J in their lattice model for river networks. 



A more realistic model than the one presented here would have coupled differential 
equations for the variables h(x,t) and the flow s. Such a coupling would capture non-local 
effects present in real basins. 
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the former for a critical reading of the manuscript. This work was supported by grants 
from INFN, NASA, NATO, NSF, ONR, EPSRC, The Fulbright Foundation, The Petroleum 
Research Fund administered by the American Chemical Society and the Center for Academic 
Computing at Penn State. AG acknowledges partial support from the HCM program under 
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FIGURES 

FIG. 1. Plot of P(s) vs. s for sizes L = 10, 30, 50, 100 evaluated at the (absolute) time t = 10 
in d = 2 + 1. The values for the parameters are: At = 0.01, Ax = 3.0, M = 1, D2/D1 = —0.1. 

FIG. 2. Collapse of the curves of P(s) from Figure 1 with r = 4/3, = (2 - t)" 1 = 3/2. 

FIG. 3. Temporal evolution of roughness G^i) = W 2 (t). The parameter used for the simula- 
tions were At = 10~ 4 , Ax = 0.5, M = 1, D 2 /Di = -0.1. 

FIG. 4. A snapshot of a typical river network created by our dynamics at two successive times 
t (a) and t' > t (b). Only sites with s > 30 are shown. The values of t and t' correspond to 89980 
and 90000 temporal iterations with At = 10~ 4 , Az = 0.5, M = 1, D 2 /Di = -0.1. The initial 
condition was a flat configuration. The flow is constructed as follows: starting from a site (x, y) 
the water can flow only to one of the three sites (x + 1, y + 1) ,(x , y + 1) or (x— 1, y + 1) thus causing 
an overall flow downward (from y to y + 1). The choice of which of the three sites is determined 
by the relative heights of the corresponding /i's (steepest descent). In the region indicated by P 
one main stream which was present in (a) has been captured by a neighboring main stream in (b) 
(piracy effect). An example of a stream present in (a) which has disappeared in (b) (abstraction 
effect) is also shown in the region A. The boundary conditions employed were periodic in the x 
direction (perpendicular to the flow) and free in the y direction (parallel to the flow). 

FIG. 5. Landscape (a) and contour plot (b) at t' corresponding to Fig.4(b). Values for the 
parameters are the same as Fig.4(b) 
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